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RESUME 

On  a mis  au  point  un  modSle  semi-empirique  de  la  propagation  de 
faisceaux  laser  en  milieux  turbulents  sous  la  forme  d'un  systSme  ferme 
d' equations  aux  deriv§es  partial les  du  second  ordre  pour  les  moments 
statistiques  du  premier  et  du  second  ordres  de  1 'amplitude  de  I'onde 
lumineuse.  Ces  Equations  sont  linSaires,  ne  contiennent  que  des  coeffi- 
cients algibriques  et  sont  uniformSment  valides  S toutes  les  distances 
de  propagation.  On  montre  que  I'hypothlse  de  la  function  de  probabilite 
normale  de  1 'amplitude  complexe  de  I'onde  est  compatible  avec  1' equation 
d'onde  stochastique.  Cette  hypothSse  est  utilisie  pour  ecrire  la  variance 
de  I'intensite  en  function  des  moments  d' ordres  un  et  deux  de  1 'amplitude. 
Les  profils  de  I'intensite  moyenne  et  de  la  variance  de  1' intensity  des 
faisceaux  laser  sont  facilement  et  rapidement  calcules  au  moyen  d'une 
methode  numerique  aux  differences  finies.  Les  resultats  sont  en  excellent 
accord  avec  des  mesures  experimental es  prises  dans  1' atmosphere  et  en 
milieu  de  simulation,  (nc] 


\ 

ABSTRACT 

A semi -empirical  model  of  laser  beam  propagation  in  turbulence 
is  written  in  the  form  of  a closed  set  of  second-order  partial  differ- 
ential equations  for  the  first-  and  second-order  amplitude  moments.  The 
equations  are  linear,  contain  only  algebraic  coefficients,  and  are  uni- 
formly valid  over  the  complete  propagation  range.  The  hypothesis  of 
normal  probability  distribution  of  the  complex  wave  amplitude  is  shown 
to  be  consistent  with  the  stochastic  wave  equation,  and  is  used  to  obtain 
the  irradiance  variance  in  terms  of  the  first-  and  second-order  amplitude 
moments.  Solutions  for  average  irradiance  and  irradiance  variance  pro- 
files of  laser  beams  are  easily  and  quickly  calculated  by  a numerical 
finite  difference  method.  Predictions  agree  very  well  with  experimental 
data  in  the  atmosphere  and  in  a simulation  medium.  (U) 
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1.0  INTRODUCTION 


Atmospheric  propagation  of  laser  beams  is  affected  by  turbulence. 
The  refractive  index  fluctuations  cause  the  beam  to  scintillate  or  break 
up  in  several  random  patches,  to  lose  its  spatial  coherence,  to  wander 
about  its  axis,  and  to  spread  out.  These  behaviours  are  detrimental  to 
military  applications  both  at  low  and  high  power  levels.  Some  adverse 
effects  are:  a reduction  of  signal -to-noise  ratio  in  detection,  ranging 
and  communication,  a loss  of  average  power  density  on  targets,  a lower- 
ing of  breakdown  threshold  owing  to  the  presence  of  intense  peaks  rela- 
tive to  average  irradiance,  etc.  Consequently,  the  modelling  of  optical 
and  infrared  wave  propagation  in  turbulence  is  important  for  the  design 
and  the  assessment  of  military  laser  systems. 


The  analytical  Rytov  method,  fully  documented  in  Refs.  1 and  2, 
gives  excellent  results  in  the  weak  scintillation  limit,  i.e.  at  small 
propagation  distance  and/or  turbulence  strength.  However,  it  is 
well-known  that  the  Rytov  method  fails  when  irradiance  fluctuations  reach 
a certain  level:  theory  predicts  an  unlimited  increase  of  irradiance 
variance  whereas  experiments  definitely  show  saturation.  Many  recent 
models  deal  with  this  problem  and  predict,  with  variable  degrees  of 
accuracy,  the  phenomenon  of  saturation.  A first  approach,  known  as  the 
renormalization  technique,  is  described  in  Ref.  2.  For  example,  de  Wolf 
(Refs.  3-6)  used  it  to  predict  saturation  of  the  irradiance  variance. 

A second  approach  consists  in  solving  a differential  equation  for  the 
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fourth-order  statistical  moment  of  the  complex  electric  field  closed  by 
local  application  of  the  method  of  small  perturbation.  This  method  ori- 
ginated in  the  USSR  and  is  reviewed  in  Ref.  7.  Examples  of  solutions 
showing  saturation  are  given  in  Refs.  8 to  11.  A third  approach  is  based 
on  the  extension  of  the  Huygens -Fresnel  principle  to  turbulent  media 
(Ref.  12).  Saturation  results  derived  from  various  methods  of  application 
of  this  principle  are  described  in  Refs.  13  to  17.  Finally,  Sancer  and 
Varvatsis  (Ref.  18)  have  proposed  an  independent  model  which  was  among 
the  first  to  show  saturation  of  the  irradiance  variance.  Although  their 
solution  was  relatively  simple  and  in  reasonable  agreement  with  data,  the 
constitutive  hypotheses  were  poorly  justified  (Ref.  19)  and  the  method 
was  unfortunately  not  revised  nor  pursued  any  further. 

It  is  not  our  purpose  here  to  discuss  the  merit  of  each  of  the 
theoretical  models  listed  in  the  preceding  paragraph.  However,  one  char- 
acteristic they  all  share  is  that  it  is  difficult  and  lengthy  to  compute 
numerical  results  from  the  proposed  solutions.  One  has  either  to  solve 
a four-fold  partial  differential  equation  or  to  evaluate  multiple  inte- 
grals at  each  spatial  point.  Moreover,  the  predicted  irradiance  variance 
is  generally  too  small  in  the  strong  scintillation  region.  Here,  we 
propose  a semi -empirical  model  in  the  form  of  linear  second-order  partial 
differential  equations  from  which  the  average  irradiance  and  the  irra- 
diance variance  profiles  can  be  easily  and  quickly  calculated  for  any 
beam  geometry. 
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This  work  was  performed  at  DREV  between  November  1976  and  October 
1977  under  PCN  21T05,  Propagation  of  Laser  Beams. 

2.0  THEORETICAL  BACKGROUND 

2.1  Equations  for  the  Instantaneous  Field 

The  basic  model  for  optical  and  infrared  wave  propagation  is 
described  in  Ref.  20.  We  shall  briefly  review  the  hypotheses  and  the 
derivation  steps. 

The  propagation  medium  is  assumed  to  have  a constant  relative 
magnetic  permeability  p and  a randomly  fluctuating  relative  dielectric 
constant  e.  Without  loss  of  generality,  we  neglect  absolution  so  that 
ep  = N^,  where  N is  the  refractive  index.  It  is  known  from  observations 
that  the  fluctuating  part  of  N is  much  smaller  than  its  average  or  unper- 
turbed value  n^  and,  therefore,  polarization  can  be  neglected.  Under 
these  conditions  the  electric  field  vector  E of  the  optical  wave  can  be 
treated  as  a scalar  E.  In  the  present  model,  E is  written  as  follows: 

E = A exp  [ik((^+z)  - iu)t]  , (1) 

where  A is  a complex  wave  amplitude;  k6,  a phase  function  to  be  determined; 
2,  the  propagation  distance;  k = n^u/c,  the  wave  number;  u,  the  angular 
frequency  of  the  assumed  monochromatic  source;  c,  the  speed  of  light  in 
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free  space;  t,  the  time;  and  i = Upon  substituting  relation  (1)  for 

E in  the  scalar  wave  equation  for  the  electric  field  derived  from 
Maxwell's  equations  (Ref.  21)  and  making  the  paraxial  approx-  'ion  valid 
for  propagation  in  atmospheric  turbulence,  we  obtain  after  separation: 

(2) 

(3) 

The  vector  V is  defined  as  follows: 

V = V*  , (4) 

and  7 is  the  gradient  operator  in  the  plane  transverse  to  the  direction 
of  propagation. 

Equation  (2)  is  the  paraxial  form  of  the  eikonal  equation.  Hence, 

V is  the  vector  angle  made  by  the  local  instantaneous  geometrical  rays 
and,  from  equation  (1),  it  follows  that  A is  the  complex  wave  amplitude 
defined  on  the  geometrical  phase  front.  The  operator  (It  * v.?)  indi- 
cates differentiation  along  a ray.  The  term  iA7.V  represents  the  geometrical 
effects  on  the  complex  amplitude  A of  the  converging  and  diverging  rays 
and  ^ 7^A  is  the  diffractional  contribution  which  tends  to  "diffuse"  A. 
Neglecting  the  latter  term  on  the  grounds  that  k is  large  constitutes 


the  geometrical  optics  approximation. 
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Since  the  refractive  index  N is  a random  function,  equations  (2) 
and  (3)  are  stochastic  equations.  There  is  no  known  general  method  for 
solving  the  random  functions  V and  A.  Even  if  that  could  be  achieved, 
it  would  yield  much  more  information  than  required  in  practice.  Here, 
the  stochastic  equations  (2)  and  (3)  are  used  to  obtain  differential 
equations  for  those  statistical  moments  of  V and  A that  are  of  interest. 


2.2  Equations  for  the  First-  and  Second-Order  Moments  of  the  Complex 
Wave  Amplitude 


We  confine  ourselves  to  the  first-  and  second-order  moments  of 
the  complex  wave  amplitude  which  lend  themselves  to  straightforward  and 
meaningful  physical  interpretations.  More  sophisticated  developments 
involving  covariances  and  higher-order  moments  are  contingent  to  the 
successful  verification  of  this  simpler  model. 

The  random  functions  are  written  as  sums  of  an  average  and  a 
fluctuating  part,  i.e.: 


N = 

n + n : 

0 

<n> 

= 0 

y 

(5) 

V = 

<V>  + V ; 

<y> 

= 0 

y 

(6) 

A = 

<A>  + a ; 

<a> 

= 0 

y 

(7) 

where  the  angular  brackets  denote  ensemble  averages.  The  equations  for 
the  moments  are  obtained  by  taking  the  ensemble  average  of  the  stochastic 
equations  for  V,  A,  AA  and  AA*  derived  from  governing  equations  (2)  and 
(3).  We  find  (Ref.  20): 
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V>.7^ 

y>.v^ 


<V>  = - iy<v.v> 


<A>  + i<A>V.<V> 


2k  r<A> 


= -y.<va>  + J<ay.v>  , 


<V>.7 


<V> 


<aa>  + <aa>7.<V>-  ^ <aa>  + ^ <7a.7a> 
-7.<vaa>  - 2<av>.7<A>  - <A><a7,y>  , 

<aa*>  + <aa*>7.<V>  - ^ 7. [<a*7a>  - <a7a*>] 
-7.<vaa*>  - <av>.7<A>*  - <a*v>.7<A> 

- 5<A>*<a7.v>  - i<A><a*7.v>  , 


where  the  superscript  * denotes  a complex  conjugate.  The  average 
irradiance  is  given  by: 

<I>  = <A><A>*  + <aa*> 


(8) 


(9) 


(10) 


(11) 


(12) 


The  physical  interpretation  of  the  mathematical  terms  in 
equations  (8)  to  (11)  is  discussed  in  Ref.  20.  We  will  simply  identify 
here  the  principal  mechanisms.  The  moments  <va>,  <vaa*>  and  <vaa>  are 
responsible  for  beam  spreading  and  <a7.v>  for  beam  breakup  or  loss  of 
coherence.  The  terms  multiplied  by  (i/k)  represent  the  explicit  diffrac- 
tionnal  effects.  Finally,  the  group  of  terms  [-<av>.7<A>*- J<A>*<a7.v>+c.c. ] 


UNCLASSIFIED 

7 

models  the  production  of  turbulent  irradiance,  <aa*>,  at  the  expense  of 
coherent  irradiance , <A><A>* . 

2.3  Integral  Solution  for  the  Fluctuating  Complex  Amplitude 


The  average  angle  <V>  is  only  weakly  dependent  on  the  variance 
of  V.  To  show  this,  we  rewrite  equation  (8): 

|—  <V>  +17  [<V>.<V>  + <v.v>]  = 0 . (8) 

To  a factor  of  order  unity,  we  have  from  small  perturbation  theory: 

<v.v>  = C^^  , (13) 

where  C ^ is  the  index  structure  constant  and  i.  , the  inner  scale  of 
n o 

14  2/3 

turbulence.  Typically,  for  strong  turbulence,  i.e.  C^^  = 2.5x10"  m”  ' , 

-3  -10 

and  = 10  m,  we  have  at  1 km  <v.v>  = 2.5x10  . By  comparison,  the 

average  ray  angle  <V>  is  of  the  order  of  where  r^  is  the  characteristic 

radius  of  the  beam  and  F,  the  focal  length  or  the  radius  of  the  initial 

-8 

phase  front  curvature.  For  r^  = 0.1  m and  F = 1 km,  we  have  <V><V>  = 10"  . 
Therefore,  for  practical  applications,  we  can  approximate: 
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~ <V>  + <V>.7<V>  = 0 

dZ  ' ^ , 


(15. 


For  an  initially  spherical  phase  front  with  radius  of  curvature  F, 
equation  (15)  has  the  solution: 


<V>  = 


(z  - F) 


(16) 


The  equation  for  the  fluctuating  amplitude,  a,  is  obtained  by 
subtracting  equation  (9)  from  equation  (3).  We  find,  using  equation  (16) 


3z  ® (z  - F) 


(17) 


where : 


g(z,r)  = - v.V<A>  - J<A>V.v 


-7.[va  - <va>]  + i[aV.v  - <a7.v>] 


(18) 


A solution  to  equation  (17)  with  boundary  conditions  a(0,r)=lim  a(z,r)“0 
is  given  by: 


a(z,r) 


CP  - ls.r'2 

2(z  - u)  (F  - u)  ~ 


^|g(u,s). 

-*  (19) 
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3.0  CLOSURE  PROBLEM 


The  system  of  equations  (9)  to  (II)  for  <A>,  <aa>  and  <aa*>  is 
mathematically  unclosed  insofar  as  it  contains  more  unknovms  than  equa- 
tions, namely  the  moments  <va>,  <yaa>,  <vaa*>  and  <a7.y>.  In  principle, 
one  could  derive  equations  for  these  moments  using  a procedure  analogous 
to  that  used  in  obtaining  equations  (9)  to  (11).  However,  one  would 
quickly  realize  that  this  technique  leads  to  equations  involving  still 
higher  order  and  unknown  moments.  This  is  the  classical  closure  problem 
which  is  always  encountered  in  the  treatment  of  statistical  phenomena 
governed  by  nonlinear  and/or  quasi -linear  stochastic  equations  as  those 
given  by  equations  (2)  and  (3).  Of  course,  there  is  no  exact  method  of 
solution  for  this  problem  since  the  complete  mathematical  model  contains 
an  infinite  number  of  equations.  Hence,  workable  models  require  closing 
the  hierarchy  of  equations  after  moments  of  a given  order.  This  can  be 
accomplished  only  through  approximations  regarding  the  higher  order 
moments  for  which  the  equations  have  been  left  out.  In  this  report,  we 
propose  four  relations  to  express  <ya>,  <vaa>,  <vaa*>  and  <a7.v>  as 
functions  of  <A>,  <aa>  and  <aa*>.  Moreover,  since  we  limit  the  analysis 
to  variances  only,  we  will  need  outside  information  to  relate  <7a.7a>  to 
<aa>. 

3.1  Statistical  Hypotheses 


We  summarize  under  this  section  the  principal  hypotheses  that 
will  be  used  to  derive  the  constitutive  relations  for  the  unknown  moments. 
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First,  we  assume  that  the  fluctuating  angle-of-arrival  v is 
statistically  homogeneous  and  isotropic  in  the  plane  transverse  to  the 
main  direction  of  propagation.  This  approximation  follows  from  the 
hypothesis  that  the  refractive  index  is  statistically  homogeneous  and 
isotropic  and  from  the  paraxial  approximation.  Indeed,  since  the  medium 
is  homogeneous  and  isotropic  and  since  the  average  ray  angle  <V>  is 
small,  the  rays  reaching  a plane  z have  traversed  nearly  equivalent  paths 
regardless  of  their  point  of  arrival  in  that  plane.  Hence,  the  covariance 
function  of  v should  depend  only  on  the  relative  positions  of  the  obser- 
vation points  in  plane  z. 


Second,  we  assume  that  the  random  amplitude,  a,  and  angle-of-arrival, 
V,  are  only  weakly  correlated.  Equation  (19)  shows  that  the  fluctuating 
amplitude  is  the  result  of  repeated  interactions  with  the  angle-of-arrival 
over  the  complete  propagation  path  between  0 and  z,  and  of  diffraction 
scattering  from  turbulent  regions  off  the  optical  axis.  Hence,  the  local 
amplitude  a(z,  r)  depends  on  many  processes  independent  of  the  local 
angle-of-arrival  v(z,  r)  and,  therefore,  the  hypothesis  of  weak  correla- 
tion appears  very  reasonable.  More  specifically,  we  assume  that; 


|<via2>|  << 


<v.v> 


<I>= 


(20a) 


<v.v> 


|<yiy2a2>|  << 


(20b) 


r 


I 

\ 


[ 


t 

1 


UNCLASSIFIED 

II 


<viV2aia2>  = '^YlY2^  <aia2> 

+ {|terms|  « '^Y'Y^  » (20c) 


etc. . . 


The  subscripts  1 and  2 refer  to  the  separation  points  of  the  covariance 
functions;  where  no  subscripts  are  used,  the  midpoint  between  1 and  2 is 
taken. 


Finally,  we  assume  that  the  covariance  functions  involving  the 
fluctuating  amplitude  are  quasi -homogeneous  and  quasi-isotropic  in  the 
transverse  plane,  i.e.: 


<a(zi,ri)  a(z2,r2)>  “ f^^(zi  ,Z2;ti ,r2)  ,^2;  lri-r2 1)  , (21a) 

<a(zi,ri)  a*(z2,r2)>  * .I2)  ’^2 ; |ri-r2 I ) , (21b) 


where  f and  f * are  functions  of  both  position  vectors  ri  and  tz,  and 

where  g and  g * are  functions  of  the  magnitude  of  the  vector  differ- 

aa  aa 

ence  only.  More  specifically,  we  set; 


f^^  ■ 1 <a(zi,ri)  a(zi,ri)>+  i <a(z2,r2)  a(z2,r2)>  , (22a) 


f^^.  ■ 1 <a(zi,ri)  a*(zi,ri)>  + i <a(Z2,r2)  a*(z2,r2)>.  (22b) 


^ . — ^ 
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Equations  (21)  and  (22)  constitute  a standard  approximation  often  made  in 
the  presence  of  average  gradients. 

3.2  Constitutive  Relations 

To  obtain  a relation  for  <aV.v>,  we  multiply  equation  (19)  by 
V.v  and  then  take  the  ensemble  average.  There  results,  inside  the  inte- 
gral operator,  an  expression  of  the  form: 

<g(u,s)  V.v(z,r)>  = -^V^.<Y(z,r)  Y(u,s)>yY^  <A(u,s)> 

- i YCu,s)>j  <A(u,s)> 

- V 7 :<v(z,r)  v(u,s)  a(u,s)> 

~s~r  - - - - 

+ J < ^7^.v(z,r)^  a(u,s)  Vg.v(u,s)>  . (23) 

From  our  second  hypothesis,  we  neglect  the  last  two  terms  of  equation  (23), 
i.e.  the  third-order  moments  involving  twice  the  angle-of-arrival  v,  as 
being  much  smaller  than  the  first  two  terms  which  are  of  the  order  of  the 
covariance  of  y*  Hence,  we  have: 


.7  <A(u,s)> 
-s 


<g(u,s)  7.v(z,r)>  = - 

- v(u,s)>^  <A(u,s)> 


(24) 


It  is  convenient  at  this  point  to  make  the  following  change  of 


variables : 
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[ 


(r,  s)  -*•  (r,  p) 


(25a) 


with: 


p « s - r 


(25b) 


from  where  it  follows: 


V -*-7-7  , 

~r  ~r  -p 


7 

-s 


7 

~P 


(25c) 

(25d) 


Also,  from  the  hypothesis  that  v is  statistically  homogeneous  and  isotropic 
in  the  transverse  plane  and  the  theory  of  isotropic  tensors  (Ref.  22),  we 
have: 


<v(z,r)  v(u,s)>  = fi(z,u;p2)  P P + f2(z,u;p^)  6 , 


(26) 


where  6 = e e + e e with  e and  e being  the  unit  vectors  in  the 
: -x  ~x  ~y  ~y  ~x  -y 

transverse  plane.  The  function  fj  and  f2  are  unspecified;  they  are  used 
here  to  illustrate  the  functional  dependence  on  variables  z,  u and  p. 

From  equations  (25)  and  (26),  we  deduce: 


7^.<v(z,r)  v(u,s)>  = 


^ ^ * 3f, 

^ 3p  ^ p 9p 


(27) 


I 

& 
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Upon  substituting  this  relation  in  equation  (24)  and  using  (19),  we  find 


for  <aV.v>: 


<a7.v>  = 


ik  r du  rr  f ik  (F-z)  ,1  1 

57/  " e ''f  [11777  tFTsr  “ J ^ JT  V sT 


+ 3fi  e-ys<A> 


-d~d’  <A(u,s)>. 


Since  the  scales  of  the  average  ainplitude  are  much  greater  than  those  of 
the  correlation  function  of  v or,  in  the  transverse  plane,  much  greater 
than  i/2(z-u)X,  the  quantities  <A(u,s)>  and  7g<A(u,s)>  can  be  taken  out  of 
the  integral  operators  in  equation  (28).  Then,  the  integrand  of  the  first 
term  on  the  right-hand  side  of  equation  (28)  has  the  form  of  a function 
of  p2  multiplied  by  the  vector  p.  Since  integration  is  over  the  complete 
plane  or  over  2Tr  radians,  this  first  term  vanishes  and  we  are  left  with: 


<a.7v>  = K(z)  <A> 


where : 


y(u,r4-e)>  . 


The  function  K(z)  is  independent  of  transverse  position  r because  of  the 
hypothesis  of  statistical  hcJmogeneity  of  v in  the  plane  normal  to  the 
propagation  axis. 
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An  analogous  procedure  is  used  to  derive  the  constitutive 
relation  for  <va>.  Multiplying  equation  (19)  by  v and  taking  the  ensemble 
average,  we  obtain,  under  the  integral  signs,  the  following  expression: 


<v(z,r)  g(u,s)>  = -<v(z,r)  v(u,s)>.  V <A(u,s)> 


- l<A(u,s)>  Vg.<Y(u,s)  v(z,r)> 


V .<v(u,s)  v(z,r)  a(u,s)> 


+ i <v(z,r)  a(u,s)  7 .v(u,s)> 


(31) 


Neglecting  the  third-order  moments  in  accordance  with  our  second  hypoth- 
esis and  using  equations  (25)  and  (26),  we  find: 


<v(z,r)  g(u,s)>  = - <v(z,r)  y(u,s)>.7  <A(u,s)> 


- i 


3fl  3f2 

P 7 — + 3fi 
3p  p 3p  ^ 


p <A(u,s)> 


Following  the  same  steps  as  before,  we  obtain: 

<va>  = - R(z)  7 <A>  , 


(32) 


(33) 


where : 


R(z) 


Jl! [2(z-u)  ■fef 


<v(z,r) .v(u,r+p)> 


(34) 
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Use  was  made  of  equation  (26)  to  reduce  the  proportionality  function  R(z) 
to  a scalar. 

The  same  derivation  technique  is  used  for  the  remaining  terms 
<vaa>  and  <vaa*>.  After  multiplication  of  equation  (19)  by  va,  the 
expression  under  the  integral  operator  becomes: 

<v(z,r)  a(z,r)  a(u,s)>  = -<v(z,r)  a(z,r)  v(u,s)>  .V^<A(u,s)> 

- i ?s."'Y(u.s)  v(z,r)  a(z,r)>  <A(u,s)> 

- V^.<v(u,s)  a(u,s)  v(z,r)  a(z,r)> 

+ i<v(z,r)  a(z,r)  a(u,s)  V^.y(u,s)>  . (35) 

From  equations  (20  b and  c)  we  obtain  to  the  order  of  <v.y>  and,  after 
rearranging : 

<v(z,r)  a(z,r)  a(u,s)>  = - <v(z,r)  v(u,s)>.7^<a(u,s)  a(z,r)> 

- i <a(u,s)  a(z,r)>  ^g.<Y(u,s)  v(z,r)>  . (36) 

Using  equations  (21),  (25),  (26)  and  (36)  and  following  the  same  steps  as 
before,  we  find: 

<vaa>  = - 1 P(z)  V <aa>  , (37) 
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where : 

(38) 

Similarly,  starting  with  l[<Y(z,r)  a*(2,r)  g(u,s)>  + <v(z,r)  a(z,r)g* (u,s)>] , 
we  obtain: 

<vaa*>  = - i Q(z)  V <aa*>  , (39) 

where : 

Q(z)  - Raalj^^^^  J'd^p  exp  pzj  <v(z.r)  .v(u.r^e)> 

* gaa*(z,u;p)  (40) 

3.3  Discussion 

We  have  derived  equations  that  relate  in  closed  form  the  higher 
order  unknown  moments  of  equations  (9)  to  (11)  to  the  lower  order  ampli- 
tude moments  for  which  a solution  is  sought.  These  relations  constitute 
the  principal  result  of  this  report  and,  for  easy  reference,  they  are 
regrouped  together  below: 


<aV.y>  = - K(z)  <A>  , (29) 
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<ya>  = 

- R(z)  V <A>  , 

(33) 

<yaa>  = 

- J P(z)  V <aa>  , 

(37) 

<vaa*>= 

-i  Q(z)  ? <aa*> 

(39) 

The  proportionality  functions  are  given  by  equations  (30) , C34) , (38)  and 
(40) . 

The  approximations  defined  by  the  inequalities  (20  a to  c)  form 
the  basis  of  our  model.  If  these  were  not  satisfied,  relations  (29),  (33), 
(37)  and  (39)  would  depend  on  cross -correlation  moments  of  y and  a,  and 
closure  of  the  system  of  equations  (9)  to  (11)  would  not  be  achieved. 

The  validity  of  this  weak  correlation  hypothesis  will  eventually  have  to 
be  confirmed  by  direct  measurements.  However,  equations  (30),  (34),  (38) 
and  (40)  show  that  the  hypothesis  is  self-consistent.  Indeed,  the  pro- 
portionality functions  vary  as  the  volume  integral  of  the  covariance  of 
y which,  for  reasonable  correlation  lengths,  should  give  results  much 

I smaller  than  the  standard  deviation  of  v. 

t 

I 

I 

A second  hypothesis  concerns  isotropy  of  angle-of-arrival  and 
quasi -isotropy  of  amplitude  fluctuations  in  the  transverse  plane.  The 

I requirements  here  are  that  the  surface  integrals  of  the  odd -order  deriv- 

1 

! 

atives  of  scalar  or  second-order  tensor  covariance  functions  be  much 
smaller  than  those  of  the  even-order  derivatives.  This  is  a less  stringent 

} condition  that  does  not  require  point  by  point  verification.  Moreover, 

I 

j since  the  major  contributions  to  the  integrals  come  from  the  neighbourhood  of 


f 
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p = 0,  the  hypothesis  of  isotropy  can  be  relaxed  to  one  of  local  isotropy, 
i.e.  one  that  applies  at  a small  separation  distance  only.  Although 
convenient,  the  approximation  could  be  abandoned  without  fundamental 
difficulties.  This  would  simply  mean  the  addition  of  one  term  in  each 
of  the  constitutive  relations:  one  term  proportional  to  7<A>  in  (29)  and 
one  term  proportional  to  <A>,  <aa>  and  <aa*>  in  (33),  (37)  and  (39)  respec- 
tively. 

The  functional  relationships  between  the  higher  order  and  lower 
order  moments  in  equations  (29) , (33) , (37)  and  (39)  are  consistent  with 
the  physical  interpretation  of  the  terms  they  model.  The  beam  breakup 
term,  <aV.v>,  is  proportional  to  the  average  amplitude  and, therefore, 
it  survives  for  plane  waves,  as  indeed  it  should.  The  terms  responsible 
i for  beam  spreading,  i.e.  <ya>,  <yaa>  and  <vaa*>,  are  all  proportional  to 

I gradients  of  the  corresponding  lower  order  mean  field  quantity  and,  as 

! 

expected,  they  vanish  for  plane  waves.  Tlie  latter  relations  are  analogous 
with  expressions  successfully  used  to  model  the  turbulent  momentum,  heat 
and  mass  transport  terms  in  turbulent  shear  flows  e.g.  (Ref.  23). 

To  determine  the  functions  K(z),  R(z),  P(z)  and  Q(z) , we  need  expressions 
for  the  covariance  function  <v(z,r)  v(u,r+p)>  and  for  the  correlation 
functions  g (z,u;p)  and  g *(z,u;p).  By  assumption  of  statistical 
homogeneity  for  v and  of  quasi -homogeneity  for  a,  the  functions  <vv>, 
g and  g * are  independent  of  beam  geometry.  Hence,  it  is  necessary 
to  solve  them  for  plane  waves  only  and,  after  substitution  in  equations 


UNCLASSIFIED 

20 


(30),  (34),  (38)  and  (40),  the  resulting  functions  K(z),  R(z),  P(z)  and 
Q(z)  are  applicable  to  any  beam  geometry.  Therefore,  the  present  approach 
conveniently  dissociates  or  separates  the  turbulence  problem  from  that 
of  geometry  as  is  the  case  with  the  Extended  Huygens -Fresnel  method  of 
Ref.  12. 

Solutions  for  the  phase  and  amplitude  correlation  functions  of 
a plane  wave  are  already  known.  However,  they  are  all  specialized  to 
displacements  in  the  transverse  plane  only.  This  is  not  sufficient  for 
the  present  purpose  and  new  developments  regarding  longitudinal  depen- 
dence are  required.  Instead  of  trying  to  solve  this  problem  exactly 
here,  we  will  proceed  indirectly  and  use  simple  known  theoretical  results 

along  with  some  empirical  data.  Comparison  with  measurements  will  show 

j 

that  this  simpler  approach  is  very  satisfactory  for  practical  applications. 

I 4.0  SEMI -EMPIRICAL  MODEL 

I 1 

4.1  Proportionality  Functions 

f 

The  approach  we  will  follow  to  determine  the  proportionality 
! function  is  indirect  or  empirical . It  is  based  on  the  requirement  that 

our  solution  for  the  irradiance  variance  Oj  should  match  the  theoretical 
[ weak  perturbation  solution  in  the  limit  of  small  propagation  distance. 


Since  the  plane  wave  solution  of  equations  (9)  to  (11)  depends  on  function 
K(z)  only,  it  is  convenient  to  start  with  this  simpler  situation. 
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Substituting  equation  (29)  for  the  beam  breakup  term  in  equations  (9)  to 
(II)  and  making  the  hypothesis  of  normal  probability  distribution  for  the 
complex  wave  amplitude,  which  will  be  discussed  later  in  this  report,  we 
obtain  the  following  plane  wave  solution  valid  in  the  limit  of  small  z: 

z 

- 4 r K(z)  dz  . (41) 

•{) 

The  corresponding  perturbation  result  is  the  Rytov  formula  for  plane 
waves  given  by  (Ref.  2): 


(42) 


Equating  (41)  to  (42)  and  taking  the  z-derivative,  we  find; 

C2 

K(z)  = 0.568  k^/^  . (43) 

o 

Strictly  speaking,  equation  (43)  is  valid  in  the  weak-scintilla- 
tion  region.  However,  since  K(z)  depends  on  phase-related  functions  only, 
as  shown  by  equation  (30) , and  since  perturbation  results  regarding  phase 
statistics  are  known  to  apply  well  into  the  strong-scintillation  region 
(Refs.  24-26),  we  simply  extend  the  validity  of  equation  (43)  to  large 
values  of  propagation  distance  without  modification.  Results  will  show 
a posteriori  that  the  approximation  is  justified. 
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For  propagation  distances  where  is  the  inner  scale 

of  turbulence  and  X the  radiation  wavelength,  the  Rytov  formula  must  be 
replaced  by  the  geometrical  optics  formula  (Ref.  2),  i.e.: 


n2  t 7/3 

0 o 


z<<£2/x 

o 


from  where  we  derive; 


K(z)  = 9.60 


C2  z2 
n 

n2£  7/3 
o o 


; z«£2/x 

o 


The  connecting  region  between  geometrical  and  diffractional 
optics  is  not  modelled  by  a simple  algebraic  relation.  For  the  present 
purpose,  we  assume  an  abrupt  passage  between  regions  of  validity  of 
equations  (43)  and  (45)  and  obtain  the  point  of  separation  as  the  distance 
z=z  at  which  expressions  (43)  and  (45)  become  equal.  We  thus  find: 


z = 0.557  £2/X 
o o 


and  the  function  K(z)  is  defined  as  follows: 


K(z)  = 


C2  z2 
n 

nH  7/3 
o o 


; "1^0 


0.568  — k7/6  z5/6  . 2 > 2 

„2  0 


(47a) 


(47b) 
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This  matching  is  very  unsophisticated  but  adequate  as  the  results  will 
show. 


There  is  no  simple  small  perturbation  solution  for  finite  beams 
from  which  R(z)  could  be  determined  in  the  same  manner  as  with  K(z).  The 
procedure  adopted  is  to  derive  an  approximate  algebraic  relation  between 
R(z)  and  the  already  known  K(z)  from  examination  of  their  respective 
functional  dependence  given  by  equations  (30)  and  (34) . These  equations 
are  of  the  same  form,  i.e.  the  same  integral  operator  applied  to  the 
covariance  functions  VpV^:<v(z,r)  v(u,r+p)>  and  <Y(z,r) .v(u,r+e)>  respec- 
tively. Approximations  for  the  latter  can  be  derived  from  the  theoretical 
expression  for  the  phase  structure  function  (Ref.  1) : 


D^(p)  = 3.44 


C2  p2 
n 

n2£  1/3 
o o 


; p < 


c2 


5/3 


D^(p)  = 2.91 


; p > I 


(48a) 


(48b) 


From  the  definition  of  vector  V,  there  follows: 


C2  z 

<v.v>  = 3.44  -2 

n2£  1/3 
o o 


(49) 


lim  V V : 

p-0  'P'P 


C2  Z 
n 

nU  V3 

o o 


<v(z,r)  v(z,r+p)> 


0.72 


(50) 
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<v(z,r) .v(z,r+c)>  _l/3 


(51) 


<v.v> 


-p-p"  YCz.r+e)> 

lim  V V :<v(z,r)  v(z,r+p)> 

p->0 


(p/%) 


-7/3 


p > I 


(52) 


Of  particular  interest  here  is  the  large  difference  between  the  decay 
rates  of  the  correlation  functions  given  by  equations  (51)  and  (52) . 

These  rates  are  to  be  compared  with  the  period  of  oscillation  of  the 
complex  exponential  in  the  integral  operator  of  equations  (30)  and  (34) . 

In  the  first  case,  i.e.  equation  (30),  the  exponential  undergoes  only  a 
few  oscillations,  except  near  u^z,  before  the  correlation  function  becomes 
negligibly  small.  Exactly  the  opposite  occurs  in  the  second  case,  i.e. 
equation  (34);  the  exponential  oscillates  very  rapidly  compared  with  the 
decay  rate  of  the  correlation  function.  Assuming  that  the  power  laws 
given  by  equations  (51)  and  (52)  are  representative  of  the  true  transverse 
variation  of  the  correlation  functions  at  nonzero  longitudinal  separation, 
we  conclude  that  the  surface  integral  of  equation  (34)  scales  roughly  as 
the  square  of  the  period  of  oscillation  of  the  complex  exponential  and 
that  of  equation  (30)  as  Therefore,  we  approximate  R(z)  as  follows: 


R(z)  = K(z) 


21 

k£2 

0 


<v(z,r) .v(z,r)> 


lim  7 7 : <v(z,r)  v(z,r+p)> 
Ip^  


(53) 
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where  t is  the  longitudinal  correlation  length  of  <v(z,r) .v(u,r+e)>. 

There  are  no  theoretical  nor  experimental  data  on  t so  all  we  can  do  here 
is  to  rely  on  a posteriori  confirmation.  We  simply  assume  that  t is  pro- 
portional to  propagation  distance  and  from  equations  (49)  , (50) , and  (53) , 
we  obtain: 

R(z)  = Y ^ K(z)  , (54) 

where  y is  an  empirical  constant  to  be  determined  by  measurements. 

The  functions  P(z)  and  Q(z)  differ  from  R(z)  by  the  presence  of 
the  correlation  functions  g and  g * under  the  integral  signs  of  equa- 
tions  (38)  and  (40).  On  an  order-of -magnitude  basis,  we  assume  that  the 
transverse  scales,  for  small  separation  distances,  of  the  correlation 
functions  g and  g * are  equal  and  can  be  approximated  by  the 
amplitude-correlation  length  for  zero  longitudinal  separation.  The  latter 
initially  increases  as  passes  through  a maximum  near  z = z^  = 

”o  and  then  decreases  proportionally  to  the  coherence 

radius  = [1.45  k^z]  of  the  Mutual  Coherence  Function  (Ref.  13); 

is  the  propagation  scale  at  which  the  coherent  amplitude  <A>  has  decreased 
by  a factor  1/e.  For  this  means  that  the  surface  integrals  in 

equations  (38)  and  (40)  depend  only  weakly  on  the  lateral  variation  of 
g and  g since  the  period  of  oscillation  of  the  complex  exponential 
3ppearing  in  equations  (38)  and  (40),  i.e.  /2 (z-uyx , is  much  smaller  than 
For  Figs.  3-5  of  Ref.  13  show  that  the  amplitude-correlation 


i 

\ 
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length,  although  decreasing  with  increasing  z,  remains  of  the  order  of 
well  into  the  strong  scintillation  region.  Hence,  we  assume  that 
P(z)  and  Q(z)  remain  practically  independent  of  the  lateral  variation  of 
g and  g ^ for  all  propagation  distances  of  interest  and,  therefore,  we 

ftcL  didi^ 

can  use  equation  (53)  to  relate  P(z)  and  Q(z)  to  K(z).  The  longitudinal 
scale  t is  not  necessarily  the  same  for  P(z),  Q(z)  and  R(z)  but,  in  the 
absence  of  better  information  and  to  minimize  the  number  of  adjustable 
empirical  constants,  we  simply  make: 

P(z)  = Q(z)  = R(z)  . (55) 

The  simple  algebraic  equations  (47) , (54)  and  (55)  constitute  an 
empirical  determination  of  the  proportionality  functions  K(z),  R(z),  P(z) 
and  Q(z).  They  are  based  on  heuristic  arguments  and  the  final  justifica- 
tion must  come  from  comparison  with  measurements.  Strictly  speaking, 
equations  (47)  and  (54)  apply  to  collimated  beams  only.  They  will  have 
to  be  revised  for  focused  beams,  especially  near  the  focal  region.  This 
will  be  done  as  soon  as  we  have  sufficient  and  varied  data  to  compare  them 
with. 

There  remains  to  model  the  moment  <Va.Va>.  From  the  hypothesis  of 
quasi -homogeneity  and  quasi -isotropy,  we  have: 

g^a  . C56) 


J 
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For  z>z^,  Yura  (Ref.  13)  has  shown  that  the  characteristic  amplitude 
coherence  scale  is  proportional  to  the  coherence  radius  p^.  Assuming  that, 
for  p not  too  large,  g has  the  same  functional  dependence  as  the  Mutual 
Coherence  Function,  i.e.  exp  [-(p/p^) , and  matching  at  P=Pjj|  this 
function  to  the  quadratic  expansion  valid  near  p=0,  we  have: 


lim 

p-K) 


-C  p-1/3  p -5/3 
m o 


(57) 


where  C is  an  empirical  constant.  The  matching  radius  p^  is  unknown  but, 
given  the  weak  dependence  of  equation  (57)  on  p^,  we  should  make  only  a 
negligible  error  by  assuming  that  it  is  constant  with  respect  to  propaga- 
tion distance.  However,  p^  should  depend  on  turbulence  condition  and 
wavelength;  the  simplest  expression  that  satisfies  these  conditions  is 
P^  Xz^.  Using  equation  (57)  with  p^  = [1.45  C^  k^z/n^]"^/?  ^nd 

2 = n^i2/iVc  12/11  k7/ii,  and  regrouping  the  constants,  we  thus  obtain: 

AO  n 


<Va.7a> 


<aa> 


(58) 


Equation  (58)  is  applicable  for  only.  However,  calculations  have 

shown  that  <Va.Va>  has  little  bearing  on  numerical  results  in  the  weak- 
scintillation  region.  For  this  reason,  equation  (58)  is  used  for  all 
values  of  z without  modification.  C is  a second  empirical  constant  and, 
since  g is  complex,  C is  also  complex  although  its  imaginary  part  is 
expected  to  be  much  smaller  than  its  real  part. 
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Finally,  from  equations  C21b)  and  {22b),  we  find: 

<a*Va>  - <aVa*>  = 0 . (59) 

4,2  Model  Equations 

The  propagation  distance  is  normalized  by  z^,  i.e.  n=z/z^,  and 
the  transverse  coordinates  by  r^.  After  substitution  of  equations  (29), 
(33)  , (37) , (39)  , (58)  and  (59)  for  the  unknown  moments , we  obtain  a 
closed  set  of  simultaneous  linear  partial  differential  equations  of 
second  order  for  <A>,  <aa>  and  <aa*>.  They  are  written  here  for 
collimated  beams: 

— <A>  + I K(ti)  <A>  - b[R(n)  + i]  v2<A>  = 0 , (60) 

9 b 

— <aa>  + i T(ri)  <aa>  - j [R(n)  + i]  v2<aa> 

= 2 bR(ri)V<A>.V<A>  + K(n)<A>2  , (61) 

<aa*>  -lb  R(n)  <aa*> 

aT\ 

= 2 b R(n)V<A>.7<A>*  + K(ti)  <A><A>*  . (62) 

The  gradient  operator  is  with  respect  to  the  normalized  coordinate  r/r^. 
The  non-dimensional  functions  K(n),  R(n)  and  T(ti)  are  given  by 


A 
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the  following  simple  algebraic  functions  of  the  normalized  propagation 
distance: 


j 0.568  ; n 1 , 

(63a) 

K(n)  = j 

[ 0.568  n-'’/®  J ^ % 

(63b) 

R(n)  = Yn  , 

(64) 

T(ti)  = Cn 

(65) 

UNCLASSIFIED 
30 

<aa*aa*>  + 2<A><aa*a*> 

+ 2<A>*<a*aa>  - <aa*>^  + 2<A><A>*<aa*> 

+ <A><A><aa>*  + <A>*<A>*<aa>  . (68) 

Equation  (68)  shows  that  Oj  depends  on  amplitude  moments  of  order  three 
and  four  which  are  not  included  in  the  present  model.  To  write  Oj  in 
terms  of  <A>,  <aa>  and  <aa*>,  we  make  the  hypothesis  of  normal  probabil- 
ity distribution  of  the  random  complex  amplitude,  a.  This  approximation 
is  based  on  the  application  of  the  central  limit  theorem  to  the  solution 
for  the  fluctuating  amplitude  obtained  from  equations  (18)  and  (19)  and 
shown  here  for  collimated  beams: 


aCz,r)  = ^ 


2,-j'  “p  [iriiT 


X |-v.V<A>  - J <A>7.v  - HaV.v  - <aV.v>]  - [v.Va  - <v.7a>]| 

(69) 


In  the  weaki-scintillation  regime,  the  linear  random  terms  v.V<A>  and 
i<A>7.v  dominate.  Neglecting  diffraction  (i.e.  k-«>) , we  have  at  small 
values  of  z: 


a(z,r) 


J'  du  { -V . 


-v.7<A>  - J<A>7.v} 


(70) 


Hence,  providing  that  the  longitudinal  correlation  length  of  v is  small 
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compared  with  propagation  distance,  the  fluctuating  amplitude  results 
from  the  sum  of  many  independent  random  processes  and  the  central  limit 
theorem  suggests  that  a(z,r)  is  normally  distributed.  As  propagation 
distance  or  turbulence  strength  increases,  the  nonlinear  random  terms 
become  dominant.  The  first  of  these  terms,  ja7.v,  gives  in  the  geomet- 
rical optics  limit: 

£n  a(2,r)  J 

Thus,  the  contribution  from  jaV.v  indicates  a log-normal  distribution. 
However,  the  effect  of  the  accompanying  nonlinear  term,  v.Va,  is  additive. 
Indeed,  this  term  represents  turbulence  scattering  from  many  different 
regions  off  the  line-of-sight;  it  is  the  multiple-path  effect  described 
by  Strohbehn  et  al.  (Ref.  27).  Hence,  this  second  nonlinear  term  suggests 
a normal  distribution  for  the  complex  amplitude  and,  as  far  as  we  can 
tell,  it  is  of  the  same  order  of  magnitude  as  jaV.v.  Finally,  equation 
(69)  shows  that,  for  2X(z-u)>£^,  the  mechanisms  leading  to  fluctuations 
of  the  complex  amplitude  and  described  above  are  all  subject  to  diffrac- 
tion scattering  from  the  turbulent  regions  off  the  optical  axis.  The 
latter  effect  is  additive  and  further  supports  Gaussian  statistics  for 
the  complex  amj>litude.  Therefore,  the  turbulence  mechanisms  of  equation 
(69)  favour  in  majority  a normal  over  a log-normal  probability  distribu- 
tion for  the  complex  wave  amplitude  A defined  by  equation  (1) . 


L 


I 


du  V.v 


(71) 
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The  assumption  of  Gaussian  statistics  for  the  complex  wave  | 

amplitude  gives:  I 

i 

I 


<aa*a*>  = <a*aa> 

= 0 

(72) 

<aa*aa*>  = 2<aa*>^ 

+ <aa><aa>* 

9 

(73) 

from  where  we  obtain: 


Oj  = <aa*>^  + <aa><aa>*  + 2<A><A>*<aa*> 

+ <A><A><aa>*  + <A>*<A>*<aa>  . (74) 

The  implications  of  the  assumption  of  normal  probability  distri- 
bution for  the  complex  wave  amplitude  are  explicitly  and  completely 
specified  by  equations  (72)  and  (73) . Whether  the  amplitude  moments  of 
order  higher  than  four  can  or  cannot  be  separated  according  to  Gaussian 
statistics  is  of  no  consequence  here.  Most  importantly,  the  solutions 
to  equations  (60)  to  (62)  are  independent  of  this  or  any  other  probability 
hypothesis.  Only  the  inference  of  the  irradiance  variance  from  these 
solutions,  as  given  by  equation  (74),  is  dependent  on  the  normal  probabil- 
ity approximation. 


UNCLASSIFIED 

33 


f 


5.0  EXPERIMENTAL  VERIFICATION 


5.1  Method  of  Solution 

Equation  (60)  for  the  average  amplitude  is  linear  and  uncoupled. 
An  analytic  solution  is  easily  obtained  in  the  case  of  an  originally 
Gaussian  beam.  It  is  given  by: 

<A>  = fCn)  exp  [-rV2r2h(T))]  , (75) 


where : 


h(n)  = 


r 

1 + ibn  + 2b J" 


R(ti)  dn 


-if  <^(n)  dn 


(76) 


(77) 


R(n)  and  /C(n)  being  algebraic  functions,  the  integrals  in  equations  (76) 
and  (77)  are  straightforward. 


Equations  (61)  and  (62)  for  <aa>  and  <aa*>  are  linear  but  coupled 
with  the  average  amplitude  through  their  right-hand  side  terms.  Closed- 
form  solutions  similar  to  equation  (75)  do  not  seem  to  exist.  However, 
general  solutions  can  be  worked  out  in  terms  of  Green's  functions.  These 
solutions  involve  multiple  integrals.  For  Gaussian  beams,  the  radial 
integration  can  be  performed  analytically  but,  even  in  this  favourable 
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case,  we  are  left  with  complicated  longitudinal  integrals  that  need  to  be 
evaluated  numerically  at  each  spatial  point,  since  the  radial  dependence 
cannot  be  factored  out.  For  this  reason,  we  find  it  more  convenient  to 
solve  the  finite  difference  version  of  equations  (60)  to  (62) , which  is 
fast  and  independent  of  beam  geometry. 

5.2  Determination  of  the  Empirical  Constants 

To  verify  our  model,  a laboratory  experiment  was  designed  which 
allows  complete  control  over  the  turbulence  parameters.  The  experiment 
is  fully  described  in  Refs.  28  and  29  and  is  illustrated  in  Fig.  1.  In 
short,  the  refractive  index  turbulences  are  produced  by  creating  an 
unstable  vertical  temperature  gradient  in  a tank  filled  with  water.  This 
is  done  simply  by  heating  the  water  at  the  bottom  of  the  tank  and  cooling 
it  at  the  top.  The  tank  is  1.5  m long,  0.6  m deep  and  0.4  m wide.  The 
propagation  path  is  increased  by  folding  the  beam  lengthwise  as  shown  in 
Fig.  1.  The  Kolmogorov  inertial  subrange  model  is  well  verified  and, 
typically,  the  index  structure  constant  C^  is  10"**  m"^/^. 

Calculations  have  shown  that  the  plane  wave  approximation,  b=0, 
is  valid  for  the  on-axis  normalized  irradiance  standard  deviation, 

6 = aj/<I>,  providing  b is  small  enough,  typically  smaller  than  0.01. 

We  use  this  to  determine  C independently  of  y since  setting  b=0  in 
equations  (60)  to  (62)  eliminates  y. 
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Equations  (60)  to  (62)  with  b=0  were  solved  for  various  values  of 
the  complex  constant  C.  The  predicted  plane  wave  3 obtained  in  this  manner 
was  compared  with  the  experimental  3 measured  on  the  axis  of  a 25  mm  diam 
collimated  Gaussian  beam.  The  constant  C that  produced  the  best  overall 
fit  is  given  by: 

C = 0.31  - 0.01  i . (78) 

The  results  are  plotted  in  Fig.  2 against  the  Rytov  formula  for  the 
log-irradiance  standard  deviation,  Two  theoretical  curves  are  shown 

corresponding  to  the  extreme  values  of  as  given  in  Table  I.  The 
agreement  is  excellent;  in  particular,  the  predicted  maximum  of  about 
1.3  and  the  subsequent  slow  decay  toward  unity  (supersaturation)  are 
very  well  corroborated. 

The  constant  y is  determined  by  a best  fit  to  the  on-axis 

[ average  irradiance  data.  The  results  are  shown  in  Fig.  3 where  the 

t 

f 

I solid  curve  represents  the  theoretical  solution  using: 

I Y = 2.8  . (79) 

t 

■ 

Without  turbulence,  the  on-axis  average  irradiance  would  remain  essen- 
tially equal  to  one  for  all  propagation  distances  at  which  the  measure- 
ments shown  in  Fig.  3 were  taken.  Hence,  the  turbulence-induced  beam 
spread  is  important  and  it  is  accurately  predicted  by  our  model . 
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FIGURE  2 - The  normalized  irradiance  standard  deviation,  8=aj/<I>, 

measured  on  the  axis  of  a 25  mm  diam  collimated  laser  beam 
plotted  versus  the  Rytov  formula  for  the  log-irradiance 
standard  deviation,  The  experimental  parameters  for 

the  corresponding  symbols  are  listed  in  Table  I.  The  solid 
curves  are  theoretical  solutions  for  0^=2  (upper  curve)  and 
nQ=3.6  (lower  curve)  obtained  with  the  method  of  this  paper. 
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FIGURE  3 - The  average  irradiance  measured  on  the  axis  of  a 25  mm  diam 
collimated  laser  beam  plotted  versus  the  Rytov  formula  for 
the  log-irradiance  standard  deviation,  The  experimental 

parameters  are  listed  in  Table  I.  The  solid  curve  is  the 
theoretical  solution  obtained  with  the  method  of  this  paper. 
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5.3  Universality  of  the  Empirical  Constants 

We  compare  here  the  theoretical  solutions  with  measurements  made 
under  various  experimental  conditions.  All  calculations  were  performed 
with  the  values  of  C and  y given  by  equations  (78)  and  (79) . 

Figure  4 shows  the  on-axis  average  irradiance  measured  in  the 
simulation  medium  with  a collimated  beam  of  radius  mm.  As  expected 

from  the  smaller  value  of  the  beam  radius,  the  resulting  beam  spread  is 
more  pronounced  than  in  Fig.  3,  which  is  very  well  predicted  by  our  model. 

Examples  of  average  irradiance  profiles  in  the  weak  and  strong 
scintillation  regions  are  plotted  in  Figs.  5 and  6 respectively.  The 
corresponding  results  for  the  irradiance  standard  deviation  are  given  in 
Figs.  7 and  8.  Data  were  obtained  from  continuous  DC  and  RMS  records  of 
the  instantaneous  irradiance  signal  derived  from  a photodetector  scanning 
the  beam  very  slowly,  about  0.02  mm/s.  In  all  cases,  agreement  between 
theory  and  measurements  is  very  good.  Experimental  errors  are  approx- 
imately 10  to  15%  of  on-axis  values.  These  relatively  large  errors  are 
due  to  the  very  sharp  irradiance  peaks  that  are  characteristic  of  optical 
scintillation.  The  crest  factor  of  the  instantaneous  signal  with  respect 
to  its  RMS  value  is  typically  of  the  order  of  10  but  it  often  reaches  a 
level  as  high  as  20  to  50.  This  requires  an  unusually  wide  dynamic  range 
for  the  instruments,  especially  for  the  RMS  meters  which  had  to  be  operated 
at  or  below  25%  of  full  scale.  Consequently,  relative  errors  of  10  to  15% 
are  to  be  expected. 


d 


k 
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The  proposed  model  must  of  course  be  applicable  on  the  atmospheric 
scale.  To  verify  this,  we  use  the  atmospheric  data  reprinted  from  Fig.  19 
of  Ref.  7 with  permission  of  IEEE,  Inc.  As  is  always  the  case  in  field 
experiments,  measurements  were  taken  for  different  values  of  the  similar- 
ity parameters  and  b.  However,  both  and  b are  smaller  than  0.01, 
as  shown  in  Table  I,  and,  thus,  can  be  made  equal  to  zero  with  no  detect- 
able errors  on  the  solution  for  the  normalized  irradiance  standard  devia- 
tion 3.  The  results  are  illustrated  in  Fig.  9.  The  agreement  between 
theory  and  data  is  as  good  as  can  be  throughout  the  propagation  range. 

This  is  to  be  contrasted  with  recent  models  where  the  predicted  3 is  too 
small  in  the  strong  scintillation  region,  e.g.  Refs.  7 and  10. 

The  universality  of  the  empirical  constants  C and  y cannot  be 
established  in  all  certainty  but  the  results  of  Figs.  2 to  9 indicate 
that  they  are  very  likely  universal.  Differences  by  a factor  of  two  in 
beam  diameter  are  very  well  handled  by  the  theory.  However,  the  strongest 
test  comes  from  the  successful  prediction  of  simulated  and  atmospheric 
data.  Dimensionally,  the  all-important  governing  parameter  differs  by 
a factor  of  about  500  and  the  propagation  distances  by  about  300.  These 
represent  experimental  conditions  that  are  certainly  as  far  apart  as  can 
be  expected  and  the  model  works  extremely  well  in  both  cases,  as  shown  by 
Figs.  2 and  9. 

i 

I 
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FIGURE  4 - The  average  irradiance  measured  on  the  axis  of  an  11  mm  diam 
collimated  laser  beam  plotted  versus  the  Rytov  formula  for 
the  log-irradiance  standard  deviation,  Oj^.  The  experimental 
parameters  are  listed  in  Table  I.  The  solid  curve  is  the 
theoretical  solution  obtained  with  the  method  of  this  paper. 
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FIGURE  5 - The  average  irradiance  profile  of  a 25  mm  diam  collimated 

laser  beam  measured  at  aj^=1.2.  The  experimental  parameters 
are  listed  in  Table  I.  The  solid  curve  is  the  theoretical 
solution  obtained  with  the  method  of  this  paper. 
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FIGURE  7 - The  irradiance  standard  deviation  profile  of  a 25  mm  diam 
collimated  laser  beam  measured  at  a_=1.2.  The  experimental 
parameters  are  listed  in  Table  I.  The  solid  curve  is  the 
theoretical  solution  obtained  with  the  method  of  this  paper. 
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FIGURE  9 - The  normalized  irradiance  standard  deviation,  6=aj/<I>, 
measured  in  the  atmosphere  and  reprinted  from  Fig.  19  of 
Ref.  7 with  permission  of  IEEE,  Inc.  The  experimental 
parameters  are  listed  in  Table  I.  The  solid  curve  is  the 
theoretical  solution  obtained  with  the  method  of  this  paper. 
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The  model  equations  predict  a saturation  level  of  unity  indepen- 


dent of  geometry.  From  equations  (60),  (61),  (63b),  (65),  and  (78),  it 
is  easy  to  show  that  the  average  amplitude  and  the  amplitude  variance 
behave  asymptotically  as: 


<A>  exp  [-.155  H(r),r/r^) 

> 

(80) 

<aa>  'V  exp  [-.055  n^]  G(n,r/r^) 

(81) 

where  H and  G are  regular  functions  of  n and  r/r^.  Equations  (80)  and 
(81)  show  that  both  <A>  and  <aa>  vanish  with  increasing  n.  Therefore, 
we  have  from  equations  (12)  and  (74) : 


lim  _ <aa*> 

ri-+<»  <I>  <aa*> 


(82) 


uniformly  and  independently  of  boundary  conditions.  This  result  is  in 


good  agreement  with  a recent  theoretical  treatment  (Ref.  16)  and  seems 


to  be  confirmed  by  almost  every  set  of  experimental  data.  In  particular, 
the  predicted  slow  decay  toward  saturation  as  indicated  by  the  small 
constant  in  the  exponential  of  equation  (81)  agrees  well  with  the  exper- 


imental trend. 


i 
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Supersaturation  has  a simple  physical  interpretation;  it  occurs 
because  the  real  and  imaginary  parts  of  the  fluctuating  complex  wave 
amplitude  are  statistically  correlated  eind  have  unequal  variances.  A 
measure  of  this  effect  is  given  by  the  magnitude  of  the  moment  <aa> 
plotted  against  Oj^  in  Fig.  10  for  the  special  case  of  a plane  wave,  b=0, 
with  ng=0.  |<aa>|/<I>  rapidly  increases  from  zero  at  to  a maximum  of 

about  0.85  at  = 3 and  then  slowly  decreases  back  to  zero  asymptotically. 

That  this  behaviour  is  connected  with  supersaturation  is  clear  from 
equation  (74)  which  shows  that  |<aa>p  contributes  directly  to  a^. 


The  following  simple  example  shows  how  correlation  between  the 
real  and  imaginary  parts  of  the  random  amplitude  can  give  rise  to  super- 
saturation of  the  normalized  irradiance  variance.  Consider  the  fluctu- 
ating functions  aj  and  32  given  by; 

aj  = sin  T + i sin  t , (83) 

32  = sin  T + i cos  t , (84) 

where  t is  a random  variable  uniformly  distributed  between  0 and  2-n . The 
correlation  factor  between  the  real  and  imaginary  parts  of  aj  is  one  and 
it  is  zero  for  32.  Assuming  for  simplicity  that  <Ai>  = <A2>  = 0 and 
then  calculating  the  normalized  irradiance  variances  corresponding  to 
aj  and  32,  we  find: 
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<aiai*  aiai*> 


<aiai*>^ 


<sin^  T> 


<sin^  T>^ 


(85) 


<a2a2*  a2a2*> 


<a2a2^ 


<sin‘*  T>  + <005**  T>  + 2<sin^T  cos^  t> 

- — 2 

[<sin^  T>  + <cos^  T>] 


=1. 


Therefore,  the  one  case  with  nonzero  correlation  between  the  real  and 
imaginary  parts  gives  a larger  normalized  irradiance  variance  than  the 
other  case  with  zero  correlation. 


The  excellent  agreement  between  theory  and  data  shown  in  Figs.  2 
and  9 confirms  the  validity  of  the  expressions  (72)  and  (73)  for  the 
third-  and  fourth-order  amplitude  moments.  These  relations  were  derived 
from  the  hypothesis  of  normal  probability  distribution  of  the  complex 
amplitude  A defined  in  equation  (1).  Hence,  measurements  justify  a 
posteriori  the  latter  hypothesis  for  the  purpose  of  relating  irradiance 
variance  to  the  first-  and  second -order  amplitude  moments,  and  this  for 
the  complete  propagation  range.  Normal  probability  distribution  for 
amplitude  does  not  mean,  as  implied  in  Ref.  27,  that  the  irradiance 
should  follow  a Rice-Nakagami  or  a Rayleigh  distribution.  Indeed,  the 
latter  distributions  assume  that  the  real  and  imaginary  parts  of  the 
fluctuating  amplitude  have  equal  variances  and  are  statistically 
uncorrelated.  Should  this  be  the  case,  the  magnitude  of  <aa>  would 
remain  zero  everywhere  in  contrast  with  the  results  of  Fig.  10.  Actually, 
this  figure  shows  that  the  asymptotic  limit  of  zero  for  l<aa>t  becomes 
applicable  only  for  crj^>20. 


(86) 


an 


FIGURE  10  - The  predicted  magnitude  of  the  second-order  amplitude  moment. 


|<aa>l,  for  a plane  wave,  b=0,  with  0^=0  plotted  versus  the 

Rytov  formula  for  the  log-irradiance  standard  deviation,  a„. 
1 ^ 
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Data  from  a recent  simulation  experiment  (Ref.  30)  show  a 
supersaturation  level  as  high  as  1.7.  Such  a high  value  of  6 cannot  be 
obtained  from  equation  (74)  in  the  near  plane-wave  experiment  of  Ref.  30. 
Therefore,  the  simple  normal  probability  hypothesis  does  not  apply  in 
this  case. 


The  laboratory  simulation  of  Ref.  30  was  conducted  at  very  large 

-4  -1/3 

values  of  C^,  namely  3 to  23.5x10  m , so  that  supersaturation  occurs 
for  /Xz<£,^.  Hence,  diffraction  can  be  neglected  in  equation  (69)  up  to 
and  including  the  supersaturation  region.  Moreover,  since  propagation 
distances  are  small,  the  effects  of  turbulence  scattering  from  inhomo- 
geneities off  the  line-of-sight  are  minimal.  Therefore,  the  "log-normal" 
contributions  to  the  fluctuating  amplitude  given  by  equation  (71)  are 
only  slightly  modified  by  turbulence  and  diffraction  scatterings  and  1 

i 

become  important  in  the  vicinity  of  z=2z^  where  the  "normal"  contributions,  | 

equation  (70) , begin  to  saturate  because  of  fading  of  the  coherent  ! 

amplitude  <A>. 


The  log-normal  expression  for  the  normalized  irradiance  variance 
is  given  by  (Ref.  31) : 


(87) 


where  1<A^>|  denotes  the  magnitude  of  the  complex  quantity  <k^>.  Equation 
(87)  would  replace  (74)  if  the  probability  distribution  were  exactly 
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log-normal.  This  expression  can  give  very  large  values  of  3 owing  to  the 
presence  of  the  vanishing  average  amplitude  <A>  in  the  denominator; 
actually,  3 would  diverge  to  infinity  if  equation  (87)  were  used  for 
2 > 2z^  (Ref.  31). 

Since  the  experimental  conditions  of  Ref.  30  imply  an  amplitude 
probability  distribution  in  the  supersaturation  region  closer  to 
log-normal  than  is  the  case  in  the  atmosphere  or  in  Refs.  28  and  29,  an 
expression  combining  (74)  and  (87)  should  be  used  instead  of  (74) . This 
would  give  larger  predicted  values  of  3 in  accordance  with  measurements. 
However,  as  propagation  distance  increases  with  respect  to  diffraction 
and  turbulence  scatterings  become  important  and,  with  them,  Gaussian 
statistics  for  the  complex  amplitude.  Hence,  even  if  3 reaches  a larger 
supersaturation  level  than  in  ordinary  experiments,  it  should  still 
saturate  at  a value  of  one,  which  is  very  well  confirmed  by  the  data  of 
Ref.  30.  Therefore,  the  qualitative  model  of  amplitude  probability 
distribution  derived  from  equation  (69)  constitutes  a plausible  explanation 
for  the  unusually  high  supersaturation  results  of  Ref.  30.  These 
probability  effects  have  no  consequences  on  the  model  equations  (60)  to 
(65).  However,  the  dependence  on  Prandt]  number  of  the  weak  perturbation 
solutions  derived  in  Ref.  30  means  comparable  changes  in  our  functions 
^^(b),  J^(b)  and  T(n).  This  aspect  will  be  studied  in  subsequent  work. 
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6.0  CONCLUSION 


A semi -empirical  model  of  laser  beam  propagation  in  turbulence 
was  developed.  Solutions  for  average  irradiance  and  irradiance  variance 
profiles  are  very  well  confirmed  by  measurements  over  a wide  range  of 
turbulence  conditions.  Precision  is  better  than  with  theories  previously 
known. 


The  model  is  in  the  form  of  a closed  set  of  simultaneous 
second-order  partial  differential  equations  for  the  amplitude  moments. 

The  equations,  which  are  linear  and  contain  only  algebraic  coefficients, 
are  uniformly  valid  over  the  complete  propagation  range.  Solutions  are 
easily  calculated  for  any  beam  geometry  by  straightforward  finite 
difference  techniques  and  they  require  very  little  computation  time. 

The  algebraic  coefficients  were  derived  from  simple  approximations 
and  contain  two  universal  empirical  constants.  A full  theoretical  base 
could  be  worked  out  subject  to  new  developments  on  the  dependence  of  the 
phase  and  amplitude  covariance  functions  on  longitudinal  separation. 
However,  the  excellent  agreement  between  predictions  and  data  shows 
that  the  approximate  and  empirical  method  used  in  this  paper  is  very 
satisfactory  for  practical  applications.  The  need  for  two  empirical 
constants  is  well  compensated  for  by  the  simplicity  of  the  model. 
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